module stable
 use cylinderGrid
 real(8) M_cyl(0:IA),j2_m(0:IA),m_j2(0:IA)
 real(8) j2ksi4(0:IA),int_chi(0:IA)
 integer recalculate
end module

!======================================

Subroutine cell_ij(II,JJ,C,rotlaw,dV,dM,ksi)
! to calculate the values in II,II-1,JJ,JJ-1 cell
 use accur
 use cylinderGrid
 Implicit none
 doubleprecision C,psi,Fpsi
 doubleprecision c1,EOS
 integer II,JJ
 real(8) dV,dM,ksi,ra,thetaa
 character*(*)rotlaw
!===================================

 ra=0.5d0*(r(II)+r(II-1))
 thetaa=0.5d0*(theta(JJ)+theta(JJ-1))
 phia=0.25d0*(phi(II,JJ)+phi(II,JJ-1)+phi(II-1,JJ)+phi(II-1,JJ-1))
 ksi=ra*dsin(thetaa)
 psi=Fpsi(ksi,rotlaw(1:len_trim(rotlaw)))
 rho=EOS(phiu*(C-phia-psi),1)/rho0
 dV=(r(II)**3-r(II-1)**3)/3.d0*(dcos(theta(JJ-1))-dcos(theta(JJ)))*4.d0*pi
 dM=dV*rho

 return
end subroutine

!============================================================
!============================================================

Subroutine m_cylinder(C,rotlaw)
! To find and then use the cylinder mass m=m(ksi)
 use accur
 use stable
 Implicit none
 doubleprecision dM,dV,C,ksi
 character*(*)rotlaw
 real(8) ksi_loc(IA),dM_loc(IA)
 integer now(IA),loc

!===============================================

 M_cyl=0.d0

 now=1
 do i=1,IA
  call cell_ij(i,1,C,rotlaw,dV,dM,ksi)
  dM_loc(i)=dM
  ksi_loc(i)=ksi
 end do

 loc=1
 do
  j=minloc(ksi_loc,mask=now.le.JMax,dim=1)

  if(ksi_loc(j).gt.r(loc))then
   loc=loc+1
   M_cyl(loc)=M_cyl(loc-1)
  end if
  M_cyl(loc)=M_cyl(loc)+dM_loc(j)
  now(j)=now(j)+1
  if(now(j).le.JMax)then
   call cell_ij(j,now(j),C,rotlaw,dV,dM,ksi)
   ksi_loc(j)=ksi
   dM_loc(j)=dM
  end if
  if(minval(now,dim=1).gt.Jmax)exit
 end do

! open(1010,file='Temp.dat')
! do i=0,IA
!  write(1010,'(10es23.15)'),r(i),M_cyl(i)
!print*,i,M_cyl(i)
! end do
! close(1010)
! print*,'Done',M_cyl(IA-2:IA)
! read*

 return
end subroutine

!=======================================================================
!***********************************************************************

subroutine find_jm(rotlaw)
! to find and fix j(m)=w*ksi^2
 use accur
 use stable
 Implicit none
 character*(*)rotlaw
 real(8) Fdpsi_drcyl,Fpsi
!----------------------------------

 j2_m(0)=0.d0

 do i=1,IA !-1
  j2_m(i)=r(i)**3*dabs(Fdpsi_drcyl(r(i),rotlaw(1:len_trim(rotlaw))))
 end do
! j2_m(IA)=j2_m(IA-1) !*(r(IA)/r(IA-1))**2
 m_j2=M_cyl/M_cyl(IA)

 return
end subroutine

!========================================================================
!************************************************************************

doubleprecision Function jfixedlaw(ksi) result(chi)
!     ksi - cylinder radius, kis < 1 inside the star
!     chi=rot. potential
 use accur
 use stable
 Implicit none
 real(8) ksi,x1,ri2,rip2
 integer ii,jj
!--------------------------------------------------

 if(recalculate.eq.1)then
  recalculate=0
  do ii=1,IA-1 ! find new values if j^2/ksi^4
   x1=M_cyl(ii)/M_cyl(IA)
   jj=minloc(x1-m_j2,mask=x1-m_j2.ge.0.d0,dim=1)-1
   j2ksi4(ii)=j2_m(jj)*(m_j2(jj+1)-x1)+j2_m(jj+1)*(x1-m_j2(jj))
   j2ksi4(ii)=j2ksi4(ii)/(m_j2(jj+1)-m_j2(jj))
  end do
  j2ksi4(IA)=j2_m(IA)
  j2ksi4(1:IA)=j2ksi4(1:IA)/r(1:IA)**4
  j2ksi4(0)=j2ksi4(1)

!open(123,file='Temp.dat')
!do i=0,IA
!write(123,'(10es23.15)'),r(i),j2ksi4(i)
!end do
!close(123)
!print*,'Davai'
!read*

  int_chi(0)=0.d0
  do ii=1,IA ! find integral 1/2*int(j^2/ksi^4 dksi^2) in fixed points
   int_chi(ii)=int_chi(ii-1)+(j2ksi4(ii-1)+j2ksi4(ii))*(r(ii)**2-r(ii-1)**2)
  end do
  int_chi=int_chi*0.25d0
 end if

!-------/for recalculate=0/---------------

 if(ksi.le.0.d0)then
  chi=0.d0
  return
 elseif(ksi.ge.1.d0)then
  chi=int_chi(IA)
  return
 end if

!chi=int_chi(IA)*ksi**2
!return

 ii=minloc(ksi-r(0:IA),mask=ksi-r(0:IA).GE.0.d0,dim=1)-1
 x1=ksi**2
 ri2=r(ii)**2
 rip2=r(ii+1)**2
 chi=int_chi(ii)+0.25d0*(j2ksi4(ii)*(2*rip2-ri2-x1)+j2ksi4(ii+1)*(x1-ri2))*(x1-ri2)/(rip2-ri2)
 return

end function

!========================================================================
!************************************************************************

Subroutine tau_plot(C,rotlaw)
! To find and plot the tau=Erot/|W_grav| distribution
 use accur
 use stable
 Implicit none
 doubleprecision P,ra,thetaa,C,psi,Fpsi,Fdpsi_drcyl,EOS,powrsmart,c1
 character*(*)rotlaw

!------------------------------------

 open(123,file='tau.dat')

 do i=0,IA
  do j=0,JMax
   ra=r(i)
   thetaa=theta(j)
   phia=phi(i,j)
   psi=Fpsi(powrsmart(ra,delta)*sin(thetaa),rotlaw(1:len_trim(rotlaw)))
   rho=EOS(phiu*(C-phia-psi),1)/rho0
   if(rho.le.0.d0)cycle
   c1=-Fdpsi_drcyl(powrsmart(ra,delta)*sin(thetaa),rotlaw(1:len_trim(rotlaw)))*powrsmart(ra,delta)*sin(thetaa)
   vphi=sqrt(max(zero,c1))
!CJ=CJ+rho*vphi*powrsmart(ra,delta)*sin(thetaa)*dV
!CI=CI+rho*(powrsmart(ra,delta)*sin(thetaa))**2*dV
   P=EOS(phiu*(C-phia-psi),4)/(rho0*phiu)
!   Pip=Pip+P*dV
!   T=T+rho*vphi**2/2d0*dV
!   W=W+rho*phia/2d0*dV
!   Eth=Eth+EOS(rho0*rho,5)/(rho0*phiu)*dV
   write(123,'(10es23.15)'),ra,thetaa,0.5d0*rho*vphi**2/(3*P+rho*vphi**2)
  end do
 end do

 close(123)

 print*,'Read tau.dat'
 read*

 return
end subroutine

!=====================================================================
!*********************************************************************

Subroutine eos_check
! To check and plot EOS
 use accur
 use stable
 Implicit none
 doubleprecision P,EOS,H
!---------------------------------------

 open(134,file='eos_check.dat')
 rho=1.d7

 do while(rho.le.1.d15)
  H=EOS(rho,3)
  P=EOS(H,4)
  write(134,'(10es23.15)'),rho,P,H
  print*,rho
!  read*
  rho=rho*1.001d0
 end do
 close(134)
 stop

end subroutine

